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ABSTRACT 



A six-level baroclinic forecast model suitable for numerical 
prognosis was devised. This model includes vertical motion due to 
large-scale orographic and eddy-stress effects. The model was 
programmed in Fortran IV language and applied to a test case based 
upon initial data fields for 1200GMT, 15 December 1967. A con- 
vergent 36-hour forecast was produced. Comparisons were made with 
the verifying analyses at 12-hourly intervals. 
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1 . 



Introduction 



One of the most successful numerical forecasting procedures is that 
of F. G. Shuman (1967) and his co-workers at the National Meteorological 
Center. This is a very sophisticated forecast scheme involving the prim- 
itive equations of motion. From the progressive improvement of NMC fore- 
cast skill-score measures, there seems to be little question that the 
NMC approach must be regarded as the direction of ultimate progress in 
the field of objective prognosis. 

One of the possible demerits to the use of the NMC system, at least 
in its present transitional stage, is the vast amount of sophisticated 
computer- time and memory required. Secondly, it employs the vertical 
( Tj- coordinate system^, where dV (troposheric (T ) is defined as 

-p- -p r 



Or -~ 



r T *-p r 



( 1 ) 



where 



p^ = pressure at terrain 

p = pressure at (Jj = 0, 1/3, 2/3, 1 

p* = pressure at the tropopause 

(assumed to be a material surface). 



In addition to the short time-step used in each NMC prediction 
iteration, a considerable amount of time is spent in interpolation from 
standard pressure levels (e.g.. Fig. 1) to the C f surfaces for the pur- 
poses of the prognosis. This is followed by a final interpolation back 
to standard pressure-surface levels. Thus, it was felt here that an at- 
tempt, using vorticity methods, to initiate a quasi -hemispheric multi-level 



^This description of the Shuman (f -model is not that of his most 
recent model (1967) which includes, in addition, some stratospheric 
layers. 
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forecast system could be of real advantage to the "small-shop" forecaster. 
In this attempt the conventional FNWC-NMC grid is used, but at least 



some of the terrain effects embodied in the Shuman scheme are retained, 
even though the vorticity equation approach is used. 

In this method one expresses the horizontal wind by means of its 
usual Helmholtz expression. 



v = ik x vyj * vx <2) 

Furthermore , if neither friction nor heat sources are to be considered, 
temporarily at least, Arakawa (1962) and others have shown that the 
following pair of equations represent an energy-conserving system over 
a closed mass of air: 

+ vx-vf*f§jL (3) 

V‘fvip= (4) 

The use of the linear balance equation (4) would have its usual effect 
of eliminating gravity waves, initially at least. If terrain-induced 
friction, in its broadest sense, is to be introduced, its effect would 
impinge upon eq. (3), effecting changes primarily upon the X field, and 
to a lesser extent upon the field, since the balance equation is 

considered to be in effect throughout the forecast period. The local 
time derivative form of (4), 

V-WiPt = V 2 §t (5) 



is to be used in deriving the omega equation in section 3. 

Because of the boundary layer model which is adopted in section 2, 
it is expected that point and/or area sources of vorticity will be 
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present within the grid; therefore, the prognostic step at any time 

requires a continuous process of lateral diffusion, within the diffusion 

limits suggested by Grimminger (1941), by means of a term of the form 

KV J to be added to the right side of (3). Diffusion also serves 

to prevent eq. (3) from developing computational instability when 

1/ is 2 ^ — 

frictional considerations are injected, provided ^ (4K) 2 + 

[Richtmyer and Morton (1967)], Here A ^ is the time increment, and 
AX is the grid-mesh size. 
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2. Inclusion of large-scale friction in the vorticity equation 

If one retains the relatively simple Arakawa vorticity equation, 
while introducing friction-induced vertical motion, the resultant 
modified equation would, in the absence of orography, assume the form: 



+ VX-tff = K V Z J 



( 6 ) 



Still considering a turbulent boundary layer over a level-earth surface, 
the second term on the right side of (6) is expressible, in terms of a 
turbulence-induced vertical motion _TL , by the identity 



-P s ui - 






(7) 



If the Cressman approach (1960) to the parameterization of the frictional 
stress at gridpoints is employed, one may write 






( 8 ) 



where is actually taken from seasonal values at gridpoints after 
Rung (1963) ♦ At the top of the frictional boundary layer (indexed as 
F=1 or 2 in Fig. 1) integration of (p) gives the eddy-stress vertical 
motion component Sif as follows: 



-TL f = - | )K* VX (pc* V jF V^) 



(9) 



which, after expansion, leads to 



Af »-$/°LC*V, F J^3rV(c*V 3F )-7Z,l (10) 



In eqs. (8) and (9) the Kung drag coefficients are applied in connection 
with the geostrophic wind'^p at the top of the planetary boundary 
layer. 
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If we revert now to a terrain-variable earth, another obvious term 



arises when orography is introduced, namely 



-fi T = £?)t = 



(ID 



This version of Il t is then combined with sif of (9) so that the 
resultant maximum terrain-induced vertical velocity to be applied at 
the boundary-layer top is 



-fl — IT j + Sl 



( 12 ) 



where -O. is to be defined as: 

-n-Lo“{V - Vi*? 'f /^QVV + £ V(c*V 5F )-7Z f ll (13) 

The term qualitatively accounts for upward 

motion on the windward side of a mountain range and the reverse effect 
on the leeward side. In the computation of this term, actual values 
of pressure at terrain have been computed from consideration of the 
Berkof sky-Bertoni station elevations (1955) employing a method by means 
of which the terrain pressure is hydrostatically computed from the 
thermal character of the isobaric layer in which the terrain is embedded 
(see Appendix A). 

For the most part, the Berkof sky-Bertoni smoothed station eleva- 
tions over North America do not exceed 9000 feet. On the Asiatic con- 
tinent, however, derived terrain pressures as low as 500 mb are en- 
countered. In this regard there is a dilemma: At what level should 

Vg f be considered representative when terrain reaches such high 
levels? If p^ lies between 850 and 700 mb a choice of planetary 
boundary layer at 700 mb would appear perfectly reasonable. If p 
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exceeds 850 mb by 100 mb or so, the planetary boundary layer could 
reasonably be taken at 850 mb. 

Some upper limit to the magnitude of the f vector used in 
was considered appropriate in view of the calibration of Kung*s sea- 
level geostrophic drag coefficients. For example, even with p^ ^ 500 

mb, it is reasonable to consider that the effective updraft component 
need not exceed that of the 700-mb geostrophic wind over the generally 
rough terrain involved. Also, for simultaneous use of both -TL-p and 
-TL £ in (12), it is desirable to arrive objectively at a mutually 
compatible gradient level. In view of these considerations, the follow- 
ing scheme was devised to delineate the boundary layer: if p^ is greater 

than 875 mb, p (the boundary layer pressure) is chosen as 850 mb; if 
r 

p is equal to or less than 875 mb, p is chosen as 700 mb. The choice 
of the planetary boundary layer which employs gridpoint wind data at 
either 850 mb or 700 mb is one of the simplifying features of the general 
frictional model introduced here. Because of the preservation of the 
orographic term, V^gp • V } with p equal either to 

(850) or to p (700), it is felt that some of the main physical 
aspects of the variable atmospheric boundary layer are being preserved 
by the model, while retaining an element of convenience in a 6-level 
forecast model. Furthermore a recent study by Graham (1968) indicates 
that Kung's drag coefficient gives consistently smaller values than those 
proposed by Cressman; so their use in (8) and (9) appears reasonable with 
these geostrophic winds at the selected tops of the atmospheric boundary 
layer. 

The significance of the specification of the top of the atmospheric 
boundary layer (A.B.L.) is that this level is conventionally associated 
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with the maximum frictional vertical motion (c.f. Haltiner et al. (1963), 
Krishnamurti (1968)] which then decays in some generally exponential 



manner on either side of this level. Thus, for p ^ p , the existence of 

F 

a f rictionally-induced vertical ve locity (p) may be postulated, in 
accordance with (13) , as 



SI (^) - F (-K) Ji L0 



(14) 



Above the A.B.L. an analytic function F^(p) for jfL (p) is chosen 
so as to imply that 



(loo) - O 1 

F* Wt) =/ j 



(15) 



and that F (p) is a decay factor which decreases realistically with 
a 

decreasing pressure. The choice of F (p) was made with a careful at- 

a 

tempt at fitting the mean vertical profiles of frictionally induced 
si (p) depicted by Haltiner et al. (1963) in the case p equal to 

r 

850 mb. The analytic function selected for F (p) was 

a 




which decreases in a closely similar manner to that shown in synoptic 

evidence by Haltiner et al. (1963), and with that alluded to by Estoque 

(1957). Values of F (p) are shown in Tables la and lb. 

a 

For elevations below the A.B.L. , there is no published data to draw 

upon. Here it was desired to apply once more the limiting boundary 

conditions similar to (15) for p = 1000 mb and p = p . It is also de- 

v 

sirable to consider a function F, (p) for multiplicative association 

b 

with il L0 with a smaller rate of variation with p/p than is the 

r 

case with F (p) . Such a function, somewhat analogous to (16) and subject 
a 
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to the above constraints as (16), was chosen by limiting the bracketed 
expression in (16) to the exponent 2: 



Fb(^) =■ [i- 



A, '°° c /w c J 



(17) 



The values of (p) are also in Tables la and lb. It is evident that 



the formulation of ^ ^0 ^ 






in (7), which will appear on the right 



side of (6), depends only upon g— . The finite-difference derivative 
with respect to pressure will in turn depend upon the scheme adopted for 
the finite-differencing of 63 in the vertical, a topic to be discussed 
in connection with the solution of the omega equation. The distinction 
between t-*3 and JTL is that the former represents the result of partition- 
ing the total vertical motion into a component free of frictional ef- 
fects, whereas n is the composite vertical motion resulting from 
large-scale friction. Note that the net frictional divergence in any 
column (which may include terrain) extending from the 1000-mb level to 
the 100-mb level is designed to be zero. Also since the lateral dif- 
fusion K. V 1 J is defined to be zero on the octagonal grid boundary, 
eqs. (4) and (6) still represent an energy-preserving system of equations 
in an adiabatic atmosphere. 
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3. Derivation and solution of the diagnostic omega equation 

In order to deduce the numerical solution of (6) for \p > 
it is clear that one requires values for ^ ^ at all prognostic 
levels* To obtain these vertical derivatives the diagnostic omega 
equation must first be solved at levels k = 2 to 6. 

The omega equation is derived from a combination of the vorticity 
equation and the thermodynamic equation. 

The vorticity equation, after Arakawa (1962), is now presented as 

+ vw?* vx- w -f (18 ) 



Equation (18) is differentiated with respect to pressure and multiplied 



by f to give the form 









,(19) 



The lateral diffusion term KVj has been omitted, since its use 
will be reserved for (18). Use of the time-differentiated linear 
balance equation in (19) leads to 

v'-lrt - vt-Wn = ( 20 ) 






The adiabatic thermodynamic equation is given as 



6 V> 



+ \y ¥ + cr(oj + n)= o (an 



where 



<T- 



**““ "" 1,111 



7^e dp' 

and 0^ is also expressible in terms of geopotential meters mb ^ 
(see Appendix D) . Application of the Laplacian operator, , to 

each term of eq. (21) gives 
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+ 7 2 '[ j ( i /' ) ^ n )]+ V L (<ro+Q£i)*0 ( 22 ) 

Then 1*t can be eliminated from eqs. (20) and (22) by subtrac- 

tion. Hence eq. (23) is derived: 

v*-«rw) f ? z =• ^ J 0/i >;) - (V, §*•) 

z (23) 

+ 7 £ • V X r - V 2 ( v x ■ v $ r \ - r - VP • V Wp t - V W. 

Now, use of the approximations, 

= Vf-S7$ yt , v£\ VXf« 0,>?) (24) 

•p 2 

together with del operator applied to (21) leads to eq. (25): 

vf-Vfp, t = (25) 



Hence, the final form of the diagnostic omega equation is 

vH<S(j) = * |^J? (V, >?)] - 7 l LKvJrJ- ^(vxfa) 



hip 

■r^\ + v 4' 



_r^hl 



(26) 



V{J(%^)+7X‘V$^4- - V V/L 



Equation (26) now has terms involving partial derivatives in both 
the vertical and horizontal, requiring finite-difference schemes for 
these derivatives. Appendix B presents the scheme used for finite-dif- 
ferencing in the vertical, and Appendix C presents the scheme used for 
finite-differencing in the horizontal. Since (jJ appears in a vertical 
derivative in eq. (3), a solution of (>J at every grid point for levels 
k equal to 2 and 6 requires placing boundary conditions on omega at 
the topmost level ( k = 7) and the bottommost level (k = 1) as well as 
at the lateral boundary of the octagonal grid (Fig. 2). These conditions 
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were that omega would be equal to zero at all of these boundaries. 
Similarly, the Jacobian and Laplacian terms of all other parameters in 
(26) are set to zero at the lateral boundary. Once the terms of (26) 
are finite-dif ferenced, we obtain a Helmholtz type of equation where 



V 1 + E A- <^Jk+ 1 . < 27 > 

Here is the finite-dif ferenced result of the right side of (26), 

and Qh .-/ is presumed known. Equation (26) was solved by a 3-dimensional 
Liebmann relaxation scheme where the ( n + 1 ) th iterate for any point 
is given by 



ur = or + a_r" 



(28) 



Y\ 

Here A. is the over-relaxation coefficient and the residual R at 



the nth iteration is 



^ ^ Hy ~ V*. “ (29) 



In (28), (29), the superscript n indicates the nth iterate result- 
ing from sequential Liebmann solution procedure. The convergence 

2 6 

criterion imposed is that when R becomes less than lO for each grid- 
point, the iteration procedure ceases. For the first iteration ( t = 0) 
all omegas of the right side of (29) were set to zero, and only the 
first 2 terms on the right side of eq. (26) were used to determine . 

When these omegas were computed for levels 2 through 6, then ^ was 
determined by use of the equation of continuity in the form 

v z x . ~ - r oo) 



2 



In units of mb sec 



-1 
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With a first determination of X at the 5 central levels ( k = 2, .• 

K 

6), the diagnostic equation could now be solved for a better Lx) and 
then a better X * This Cd * X iterative procedure is repeated until 
for all gridpoints, at the 5 central levels the omegas reach the con- 
vergent value. 

For times other than t = 0, the full diagnostic equation is used 
even in the first iteration. The first guess for "X is that derived 
from the previous time step. 

The above procedures apply identically for levels k = 2 ,..., 6 . 

For the 1000-mb level a somewhat different procedure is necessary to 

obtain — — for use in connection with the continuity equation (30) 

Op* 

This procedure is outlined in Appendix B. This same procedure also 
applies when obtaining for the prognostic step at 1000 mb. 

bps 
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4. Solution of the prognostic equation 

With the vertical derivatives of omega and the divergent part of 
the wind now available for levels k = 1 through 6, the vorticity eq. (18) 

is now ready for solution for these same levels. With all the advective 

0 

and divergent terms of this equation known at time t, and after finite- 
differencing and transposition, eq. (18) leads to the Poisson-type 
differential equation 

— G-(hfj-b) 



where G (i, j, k) represents the three-dimensional forcing function. 
This equation is solvable one level at a time by assuming both that 
at the octagonal boundary is zero and that all terms comprising G (i, 
j, k) are also individually zero on the boundary. For solution of (31) 
one proceeds to a Liebmann two-dimensional iterative relaxation pro- 
cedure, using successive iterates in the form 

_ , At* * -nr\ 

+• AR (32) 

A 



for any one level. Here, as in the relaxation process for (jJ , the 
superscript n + 1 refers to the (n + l)th iterative value, and the n- 
superscript refers to the nth iterative value. Again A refers to the 
over-relaxation coefficient and the residual R is expressed as: 






This iterative procedure continues until the maximum residual at any 

iteration is less than the prescribed epsilon value, € = 6.4 x 10 

-2 

sec 

Because of the simplified method of determining the gradient level, 
and to take into account discontinuities in the A.B.L. occurring between 
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adjacent gridpoints, the term K7*7 is included in (18) as a 

smoothing, or diffusion term. Several values of K were used with vary- 

5 2-1 

ing results. The initial value was 4.0 x 10 m sec , as suggested by 

Danard (1966). This resulted in extreme deepening of lows and building 

% 

6 2 — 1 

of highs. The value 4.0 x 10 m sec was then attempted which gave 

the opposite effect, namely to flatten out systems excessively. Also 

lateral boundary difficulties became apparent with the larger smoothing 

6 2 — 1 

coefficient. The next and final value chosen was 1.2 x 10 m sec which, 
while perhaps not the optimum value, gave better verification at t = 12 
hours. In accord with convergence requirements for a parabolic-type 
differential equation, the smoothing was applied to the field of 

the time t-1 rather than time t. 

With ft ( t = 0) computed initially, for a 30 minute time- 

increment was determined by a forward time step, represented symbolical- 



ly by 



= V'+ft (At) 



(34) 



For the next time step, also for ^ t = 30 minutes, a centered or 
leap-frog type time step was utilized with m = 1 in (35): 
wi+l Nll m-/ ^ m 



Y""' + WUt t) 



(35) 



For m 1 all time steps were of the leap-frog type with 2 ^ t = 

2 hours. 

New values for the stream function were thus obtained at all levels 
up to 200 mb. These values were converted to heights by the solution of 
the balance equation. In order to obtain vertical derivatives of geo- 
potential at 200 mb, and to obtain a more accurate stability parameter 



/ 
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(f (200), it was deemed necessary to obtain a good estimate of the 
100-mb height field (see (fj^ in eq. D-6). This was accomplished by use 
of regression equations developed by Lea (1961), for details see Ap- 
pendix E. 

With forecast values of Z available for levels k = 1 to 7, further 
time steps could now be computed by the procedure described in the 
diagnostic section and this section. 
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5. Forecast verification results 

As a preface to the discussion of the forecast verification results, 
the labeling of the maps (figs. 5-30) will be explained. All 1000-mb 
contours are labeled in whole meters and contours are drawn at every 
60 meters. The 500- and 300-mb contours are labeled with the thousands 
digit dropped. (Thus a height of 5580 meters is labeled 580). The con- 
tours at these levels are drawn at every 120 meters. The isolines of 

-3 

vertical velocity are drawn with an interval of 2 x 10 millibars per 

second with the exception of the 200-mb level where contours are drawn 

-3 3 

at every 0.5 x 10 interval. Labels are millibars per second x 10 . 

Although forecast maps were produced for 7 levels, only the 1000- , 
500- , and 300-mb levels are presented in order to maintain the discus- 
sion at a reasonable level. These forecasts are for 12-hour intervals 
out to 36 hours. 

line test case chosen, 1200Z 15 Dec. 1967, was one which was recent 
and thus readily available and one which shows extensive deepening and 
movement of pressure systems at most levels. 

The initial 1000-mb analysis shows an inactive 60-meter low center 
over Southern and Baja California which in 36 hours deepened to a minus 
120-meter contour and moved to the Colorado area. The 36-hour forecast 
did not verify the deepening; however, a trough-type intrusion into this 
area has occurred in the same general area. This low system aloft had 
a near vertical closed circulation throughout the atmosphere, and had a 
movement to the northeast of 600 nautical miles, with slight deepening 
at the 300-mb level. The upper level forecast moved this system north- 
ward 150 nautical miles and filled it considerably. 
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In the North Atlantic a low system of minus 120-meters, with 
troughing to the WNW at 1000 mb at initial time, split into two minus 
120-meter systems in 36 hours: one positioned near the original low 

and the second positioned just west of Newfoundland. The 36-hour fore- 
cast for this level shows a minus 240-meter contour just to the north- 
west of Newfoundland. The initial system had a WNW tilt with height 
and a closed 8640-meter contour centered at 52N-65W at the 300-mb level, 
which in 36 hours deepened to a 8520-meter contour at 47N-60W. For the 
forecast system at 300-mb the model verified on position but filled the 
low center to a 8700-meter contour. At the 500-mb level the prognosis 
verified with respect to position and height. 

Initially, east of Hokkaido, Japan, there was a 1000-mb low center 
with a minus 60-meter contour. This low deepened rapidly for the first 
12 hours and then the deepening subsided. In 36 hours it was positioned 
east of the Kamchatka Peninsula with a closed minus 180-meter contour. 

The surface forecast position verified exactly, but the forecast deepened 
the system too much. The initial 300-mb map indicated that a short wave 
trough was associated with this surface system; in 36 hours it had moved 
NE an amount consistent with that of the surface system while leaving a 
stationary long wave secondary trough west of Kamchatka. Both of these 
300-mb trough features verified quite well. 

At the initial time, a surface low center (minus 120-meter con- 
tour), was located off the central coast of Norway. It traveled ESE 
to Western Russia, maintaining its initial contour height. The model 
captured this movement almost exactly but again indicated too much low- 
level deepening. At the upper levels a trough associated with this 
system moved ESE to a position over the Baltic Sea in 36 hours. This 
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trough movement was captured well by the model at both the 500- and 300- 
mb levels and the height of the 500-mb level also verified. 

Centered at 67N-85E at the 1000-mb level at the initial time was a 
minus 180-meter contour. This low filled to a minus 60-meter contour 
in 36 hours and was positioned at 74N-110E in 36 hours. The model missed 
this system badly by deepening it throughout the 36— hour period to a 
minus 240 meter contour positioned at 72N-90E, about 250 miles to the 
west of the actual position. This system at 500 mb was only forecast 
to deepen slightly. This slight deepening and position was verified by 
the 00Z 17 December analysis. 

After 36 hours, with the exception of the low system over the 
central United States, there was generally good agreement in the 1000-mb 
low systems* positions; however, the model tested produced excessive 
deepening in some of the surface low systems by an extent of the order 
of 15 mb. At the upper levels one may easily identify comparable systems 
but again the position of systems verified better than the intensities; 
and the 500-mb systems verified better than those of the 300-mb level on 
both intensity and position. Table 2 presents the pillow and root mean 
square error for all 6 levels at 12 -hourly intervals based on all in- 
ternal gridpoints . 

Figures 26 to 30 show the vertical-motion or Gj charts produced 
at t = 12 hours, namely at 00Z, 16 December 1967 for the levels 850, 

700, 500, 300, and 200 mb. It is noteworthy that maximum values seem 
to be closely aligned in the vertical with extreme maximum values oc- 
curring generally at 500 mb. However, U 300 was in general closer in 

magnitude to ^-^qq than was The extreme magnitudes computed were 

-3 -1 

near 6.0 x 10 mb sec at 500 mb, which are close to values reported 
in a case study by Krishnamurti (1968). 
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6 . 



Remarks and conclusions 



Due to lack of sufficient core space it was not possible to main- 
tain all necessary fields for one time step in core. The choice was made 
to maintain only the information for a triad of levels at one time. This 
necessitated a large amount of reading from and writing on direct-access 
storage devices as is evidenced by the FORTRAN program. This increased 
the running time of the program considerably. In order to limit the 
continuous running time of the program, the 36-hour forecast was pro- 
duced in 3 runs of 12 hours each. The current values of height, 6 l) , X > 

as well as the at the preceding time step were stored on disk 

and these values were read to begin the first time step of the next run. 
(For the layout of direct-access devices see figures 4a, 4b, and 4c.). 
Running time was 110 minutes for the first 12 hours and 90 minutes for 
each of the remaining 2 runs. This time could be reduced to a more 
operationally acceptable figure by changing from Fortran to machine 
language programming, conversion to one-dimensional arrays, and, in 
general, increased programming efficiency. 

As originally planned, the 100-mb level was to have been completely 
inactive; that is, the equation 

+ VyrVYl + VX* VP -f ^ + S\~) (36) 

as applied at 100 mb is subject to the following boundary conditions: 

%= o 

Vy *V^ = O 
X = 0 
6J+ .a = 0 
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In order to deduce improved values of vertical derivatives of the 



geopotential at 200 mb thus improving the static stability parameter, 
the heights at the 100-mb level were altered at each time step using the 
Navy regression procedure outlined in Appendix E. This procedure was in- 
tended to provide a reasonable estimate of the mean temperature in the 
200 to 100-mb layer for use in the thermodynamic equation* Due to an 
inadvertent programming error, however, the (§) values were converted 
to a *l// ^oo which varied in time. This lapse caused the altera- 

tion of the static boundary conditions above and, through the equa- 
tion, it affected the values of -p in eq. (26), especial- 

ly at the upper levels. The actual desired result of this term at 200 
mb in centered finite-difference form is 



(37) 

200 Yn~b“ 

The effects on this important term of the CJ - equation were opera- 
tive from t - 0 but were not discernible until diverging values of 
were detected at t = 42 hours. It is planned to rectify this problem and 
report the results at a future date. 

As is depicted in Table 2, the model appears to produce too much 
deepening at low levels and an overall increase in heights above 850 mb. 
The reason for this phenomenon is as yet undetermined. The fact that 
systems do indeed maintain their identity and that flow patterns, at least 
in a qualitative sense, show much similarity to the verifying analyses 
indicates that the model has some merit and further work on its refine- 
ment is justified. 
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TABLE 1 



F (p) VALUES FOR USE IN EQUATION 14 
(a) Gradient Level Assumed as 850 mb 



Pressure level 


F(p) 


1000 mb 


0.00000 


850 mb 


1.00000 


700 mb 


0.75178 


500 mb 


0.42534 


300 mb 


0.13528 


200 mb 


0.03398 


100 mb 


0.00000 


(b) Gradient Level Assumed as 700 mb 


Pressure level 


F(p) 


1000 mb 


0.00000 


850 mb 


0.20763 


700 mb 


1.00000 


500 mb 


0.56579 


300 mb 


0.17995 


200 mb 


0.04520 


100 mb 


0.00000 
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Level Diagram 




k = 7 
k = 6 

k = 5 



k = 4 



k = 3 



k = 2 



k = 1 



F 



Schematic representation of the pressure surfaces 
used in the diagnostic procedures. The index 
k = 1,...,7 is identified with the constant pressure 
values shown in mb. The index F - 1, 2 denotes the 
top of the atmospheric boundary level or gradient 
level, of 850 mb, or 700 mb, the selection of which 
depends upon the station pressure p r ^ (i,j ). 
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Figure 2 
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FLOW DIAGRAM 





Figure 3a. 



FLOW DIAGRAM (continued) 




Figure 3b 
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TAPE 11 RECORD LAYOUT 



Record 

1 . 

2 . 

3. 

4. 

5. 

6 . 

7. 

8 . 

9. 



1000-mb initial 

850-mb initial 

700-mb initial 

500-mb initial 

300-rab initial 

200-mb initial 

100-mb initial 

2 -values 
o 

Z^-values 

Figure < 



Contents 
Z values 
Z values 
Z values 
Z values 
Z values 
Z values 
Z values 
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DISK 12 AND 13 RECORD LAYOUT 



Record 

1 . 

2 . 

3. 

4 . 

5. 

6 - 10 . 

11-15. 

16-20. 

21-25. 

26-30. 

31-35. 



Contents 

1000-mb Z intermediate forecast 
1000-mb previous time intermediate forecast 

1000-mb 'p current time intermediate forecast 

1000-mb (0 intermediate forecast 

1000-mb intermediate forecast 



contents 


& 


sequence 


same 


as 


Records 


1-5 


for 


850 


mb 


contents 


& 


sequence 


same 


as 


Records 


1-5 


for 


700 


mb 


contents 


& 


sequence 


same 


as 


Records 


1-5 


for 


500 


mb 


contents 


& 


sequence 


same 


as 


Records 


1-5 


for 


300 


mb 


contents 


6c 


sequence 


same 


as 


Records 


1-5 


for 


200 


mb 


contents 


6c 


sequence 


same 


as 


Records 


1-5 


for 


100 


mb 



Figure 4b 
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TAPE 14 RECORD LAYOUT 



Record Contents 

1. 1000-mb Z final forecast"^ 

2. 1000-mb previous time final forecast 

3. 1000-mb current time final forecast 

4* 1000-mb (j) final forecast 

5. 1000-mb final forecast 



6-10. 


Contents 


6c 


sequence 


same 


as 


Records 


1-5 


for 


850 


mb 


11-15. 


Contents 


6c 


sequence 


same 


as 


Records 


1-5 


for 


700 


mb 


16-20. 


Contents 


6c 


sequence 


same 


as 


Records 


1-5 


for 


500 


mb 


21-25. 


Contents 


6c 


sequence 


same 


as 


Records 


1-5 


for 


300 


mb 


26-30. 


Contents 


6c 


sequence 


same 


as 


Records 


1-5 


for 


200 


mb 


31-35. 


Contents 


6c 


sequence 


same 


as 


Records 


1-5 


for 


100 


mb 


36-70. 


Contents 


6c 


sequence 


same 


as 


Records 


1-35 


for 2nd 



forecast interval 



Figure 4c 



o 

’’Note : Although these records are designated as final forecast, 

they are used as input for the next computer interval run if the total 
forecast is partitioned into forecast interval segments. 
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Figure 5 
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Figure 6 
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Figure 7 
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Figure 8 
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Figure 9 
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Figure 10 
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Figure 12 
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Figure 13 
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Figure 14 
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Figure 15 
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Figure 16 



53 




Figure 17 



54 - 




Figure 18 



-55 




56 - 







/ 



r ~ 



r 

<k 



t 



Figure 20 



- 57 - 




• 58 - 




Figure 22 
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Figure 23 
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Figure 24 
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Figure 25 
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Figure 26 
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Figure 27 
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Figure 28 
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Figure 29 
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Figure 30 
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APPENDIX A 



Computation of Terrain Pressure 



In exhibiting the method of determination of terrain pressure 
the example used will be the case of terrain lower than the 850-mb 
standard height* Here is the terrain height, is the 1000-mb 
height and Z^ is the 850-mb height, We start with the hypsometric 
relationship: 



P 



T 




8 



T JsVxJ 



(A-l) 



In (A-l), T is the mean temperature in the layer. Equation A-l may also 
be expressed as 



A2 = A Z iT D+ Sj 



(A-2) 



where is the mean specific virtual temperature anomaly as dis- 

cussed in Haltiner and Martin (1957) and subscript ST refers to stand- 
ard atmosphere value. After replacing temperature. 



T- _ J_ d $ 

"c) JL™ — — 

the following approximation for 1 + Sv is obtained in finite-differ- 



ence form: 



1 +- Sv = $ h Z 1 1 -- 

($z - $>) ST 






(A-3) 



Since A Z is known at any time, x (2 T ~^')sT ma ^ ^ eter “ 

mined by the substitution of (A-3) into (A-2); hence, a value of 
may be obtained. This value is compared with standard atmosphere heights 
from U. S. Standard Atmosphere (1962). When the height values are matched, 
the corrected terrain pressure is determined. 
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For cases in which the terrain extends above the 850-mb standard 
height, the same procedure is followed and sub 1 refers to 850 mb and 
sub 2 refers to 700 mb, and similarly for higher terrain heights. 
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APPENDIX B 



Treatment of Vertical Derivatives 

For fitting tO (vertical motion) profiles to three successive 
pressure levels (lower level subscript 1, middle level subscript 2, 
upper level subscript 3) a Lagrangian quadratic interpolation formula 
is used on omega. 

v(p)~ )(r-h) + {f- rjifcrsk f£j3iS£^t . i ) 

After differentiating the above and setting p = the following equa- 
tion is obtained: 

^ -■ >, m-M- + , ( b - 2 ) 

«-*)( im §) w-ny-ii-ty 

Similarly s differentiating (B-l) twice with respect to pressure and set- 
ting p = p^ the following equation is obtained: 

] — pf _ i ^± -r H ^.3.-- *~\(B-3) 

~ L cn-tXi-ty (^xi-ny 

sl 9 the frictional and terrain induced vertical motion, is treated 
the same as £J in relation to the format of vertical derivatives. 

For fitting ^ (geopotential) and \fJ (streamfunction) to three 
successive pressure levels the same procedure, as above is followed with 
the exception that (^ and \j) are considered to have a quadratic 
vertical variation on a logarithmic pressure scale leading to the follow 
ing results for centered at p = p^. 
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$(?)= 1 , 






+ 



+ * /P-TOft-Pai 

HV>i)(VP 5 ) ^ ft-'RjCPj-TV) <B ’ 4) 



where P = Ln(p) and 



3$ \ _ j_ riAi^L- + hJk (. &) ~u :) 

F° r at t * ie lOOO-mb level, omega is expressed as a cubic equation 

in terms of omega at the 1000-mb (which is equal to zero), 850-mb, 700- 
mb, and 500-mb levels. After differentiating with respect to pressure 
and setting p = P-^qqq the following equation is obtained: 



d-p 



ps.|0oo 



(?<-*>) Lt-fyCfi-Pi) 

&--PJ **R)+f. A 

ti 

<&-T.)Lh-*)(.*3-*L) 



^?So + 
^ 7.0 t (b- 6 ) 
o 



Setting in appropriate values of P^, P^, P^, and P , the above becomes, 

b^\ - olio H 1C, . oog35S3 ^ 2 . ( B - ? ) 

Z>p)r»--< “» -.001X857^3 

The coefficients derived above are constant with respect to time. 
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The values of the coefficients appearing in 



$or) 



and 






given by eqs* ( 3 - 4)5 (B-5) have been stored and employed in subroutine 
for generating 



$ ( 1 k - l 1 

\ k - 2 6 

£> *> 1 % 

In the uppermost level, i.e. for both $ 7 and ( )g it was found 

to be useful to make use of the Navy Weather Research Facility extra- 
polation procedure to obtain a hydrostatically-consistent best estimate 

£or $ 100- 
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APPENDIX C 



Treatment of Horizontal Derivatives 

The horizontal derivatives this paper is concerned with are the 
Laplacian operator, the Jacobian operator and the del-dot-del opera- 
tion. The finite difference forms of these are exhibited below; the 
superscripts refer to the designated points exhibited in Fig. 31. 

Laplacian: v ' + + ^ Ae) (C-l) 

Jacobian: J>/S) [( B,)- ( A: (8,- Bj)J (c- 2 ) 

Del A-Del B - W7B = ^ [(>V (fif ■ (&- B^)j «=-3 ) 

The map factor 'W* in the above equations is dependent on latitude (fa 
and is expressed as: 



y^s i» sin at 

1 + S \Y\ 60 



(C-4) 



The distance between gridpoints, "d" in the grid used (Figs. 2, 
31), is 381 kilometers, true at 60 north. 
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4 



Finite Differencing Grid 
Figure 31 
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APPENDIX D 



Computations of Static Stability Parameter 



The static stability parameter, CT , at level k may be expressed 
as follows: 




LliZL aa') 

\ 0 d 



which may be rewritten as: 



4 = 

by use of the Poisson relationship. Using the relation 



(D-l) 



(D-2) 



T=. _J- 

gives: 

d a i./ii _ Jk M \ 

** r Vdf.M*') 1 C ^ dSU-p) 



(D-3) 



(D-4) 



Assuming constant within the layer pair (k-1, k) , (k, k+1) 

the previous equation is set up as an ordinary differential equation in 

the form 



ot^ _ 7U dj. 

C X Cfcnl*)* C -p d Jtn'P 



4?^ 



(D-5) 



The above is integrated with respect to JlrJ p from ^ 
where these pressures are the log means for the two layers. 



e ° p kH 

that is: 
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After the above operation has been performed and solved for (fj^ , it 
leads to the finite-difference form: 







JuH± zc*, 

A 



(D-6) 



The lower bound on the braced quantity above is to be 
thus ma Y solved at all required levels* 
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APPENDIX E 



Procedure for 100-mb Height Extrapolation 

Navy Weather Research Facility (Lea, 1961) outlines a method of 
deriving a 100-mb extrapolated height field from the 200-mb heights 
and temperature fields in the form 



A 0 + ^\^-~2oo* ^Z~^2oo 



(E-l) 



The coefficients A^, A , and A^ vary with both latitude and month of 
the year. 

The regression coefficients Aq, A^, A^ are defined for the centers 

o 

of latitude belts of 10 latitude (or multiple thereof) in width. Since 
the A^ for one belt generally is different from that of the adjacent 
belt, continuity required a smooth transition from one belt into the 
next, which was accomplished by linear interpolation between A^' s and 
the linear distance between latitude band centers. 

Since temperature fields were not utilized directly in this study, 
the following approximation was made! 



j_ 

1U ZJUvp' 



(E-2) 



Using the centered finite-difference scheme for 
equation (E-l) becomes 



dj. 

bln-p 



at 200 mb. 



/ $3*Q 

Z — A A t A 2 . v’R e ^.£vi3 

IOO 



273) 



(E-3) 



(l-M* ~) 



79 



FORTRAN IV PROGRAM 
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